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Abstract. MHD turbulence plays a central role in the physics of star-forming molec- 
ular clouds and the interstellar medium. MHD turbulence in molecular clouds must be 
driven to account for the observed supersonic motions in the clouds, as even strongly 
magnetized turbulence decays quickly. Driven MHD turbulence can globally support 
gravitationally unstable regions, but local collapse inevitably occurs. Differences in the 
strength of driving and the gas density may explain the very different rates of star 
formation observed in different galaxies. Two types of comparisons to observations 
are reviewed. First, the use of wavelet transform methods suggest that the driving 
comes from scales larger than observed molecular clouds. Second, comparison of simu- 
lated spectral cubes from models to real observations suggests that Larson's mass-size 
relationship is an observational artifact. The driving mechanism for the turbulence 
is likely a combination of field supernovae in star-forming sections of galactic disks, 
and magnetorotational instabilities in outer disks and low surface brightness galaxies. 
Supernova-driven turbulence has a broad range of pressures with a roughly log-normal 
distribution. High-pressure, cold regions can be formed even in the absence of self- 
gravity. 



1 Introduction 

One of the big questions in star formation is what determines the rate of star 
formation in galaxies? Another, more pointed way of phrasing this question is 
to ask why the star formation rate in normal galaxies is so low, and why it 
varies so strongly, over orders of magnitude from low surface brightness galax- 
ies, through normal galaxies, to starburst galaxies. The free-fall time for gas at 
typical interstellar densities is 

where p is the mean mass density of the gas, G the gravitational constant and 
n = the number density, with p = 2.36mH- Yet galactic ages range up to 
10^° yr, and star formation continues today. What has delayed star formation 
sufficiently to allow it to continue? 

In what might be called the standard theory of star formation, magnetic 
fields are invoked to answer both of these questions. If fields are strong enough, 
they can magnetostatically support clouds against collapse. The star formation 
rate would then be determined by the rate of ambipolar drift of neutral gas 
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past ions tied to the magnetic field towards the centers of self-gravitating cores 
[^2|7^. Furthermore, if the fields are strong enough that the Alfvcn speed va 
reaches the rms velocity v, then strong shocks will be converted to MHD waves. 
As linear Alfven waves are lossless, it was thought that motions remaining from 
the initial formation of the clouds might be enough to explain the observation of 
strongly supersonic motions in molecular clouds Q . In this review I will explain 
why both of these ideas now appear questionable. 



2 Decaying Turbulence 

First let us consider the decay of supersonic turbulence as shown in Figure |l[ 
I call on computations performed with two different methods: Eulerian hydro- 
dynamics and MHD on a grid, using the code ZEUS-3D |]83| , p4|j2^ ,|3^ , available 
from the Laboratory for Computational Astrophysics at 



http : / /zeus . ncsa. uiuc . edu/lca-home_pagc . html 



and Lagrangian hydrodynamics using a smoothed particle hydrodynamics (SPH) 



code derived from that described by Benz |21 and Monaghan running on 
special purpose GRAPE processors |^,|2| , and incorporating sink particles . 

We chose initial conditions for our models inspired by the popular idea that 
setting up velocity perturbations with an initial power spectrum P{k) cx in 
Fourier space similar to that of developed turbulence would be in some way 
equivalent to starting with developed turbulence j6|,^,Q. Observing the de- 
velopment of our models, it became clear to us that the loss of phase information 
in the power spectrum |^ allows extremely different gas distributions to have 
the same power spectrum. This is particularly important for supersonic flows. 
Supersonic, HD turbulence has been found in simulations [ fz^ to have a power 
spectrum a = —2. However, any single, discontinuous shock wave will also have 
such a power spectrum, as that is simply the Fourier transform of a step func- 
tion, and taking the Fourier transform of many shocks will not change this power 
law. Nevertheless, most distributions with a ~ —2 do not contain shocks. 

After experimentation, we decided that the quickest way to generate fully de- 
veloped turbulence was with perturbations having a flat power spectrum a = 
for < fee; < 8. We set up velocity perturbutions drawn from a Gaussian random 
field fully determined by its power spectrum in Fourier space following the stan- 
dard procedure: for each wavenumber we randomly select an amplitude from 
a Gaussian distribution centered on zero and with width P{kd) = Pok^ with 
kd = |kd|, and a phase between zero and 27r. We then transform the field back 
into real space to obtain the velocity in each zone. This is done independently 
for each velocity component. For the SPH calculation the velocities defined on 
the grid are assigned onto individual particles using the "cloud-in-cell" scheme 
[ p] . In all of our models we take Cs — 0.1, initial density po — 1, and we use a 
periodic grid with sides L = 2 centered on the origin. These parameter choices 
define our unit system. Our choice of periodic boundary conditions corresponds 
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Fig. 1. Log of density is shown at times t/ts — 0.2, 0.6 and 1.0 on slices through 
the decaying supersonic hydrodynamic models C and D with initial rms Mach number 
M = 5 described in jsHj at standard resolution (128"^) and high resolution (256^), where 
ts is the sound crossing time of our numerical box. Note that each image is scaled to 
its own maximum and minimum to enhance morphological features (From p^) 

to the case of free turbulence discussed above, at least initially. Thereafter, the 
appropriate treatment is less clear. 

We next performed resolution studies using ZEUS for three different cases 
with no field, weak field and strong field as shown in Fig. |^. The weak field 
models have an initial ratio of thermal to magnetic pressure P — 2, while the 
strong field models have f3 = 0.08. We ran the same unmagnetized model with 
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Fig. 2. Three-dimensional resolution studies for initial rms Mach number M = 5, 
isothermal models. ZEUS models have 32'' (dotted), 64"' [short dashed), 128^ (long 
dashed), or 256'^ (solid) zones, while the SPH models have 7000 [dotted), 50,000 (short 
dashed), or 350,000 (solid) particles. Panels show a) hydro runs with ZEUS, b) hydro 
runs with SPH, c) initial rms Alfven number A = v/va = 5 MHD runs with ZEUS, 
a.nd d) A=l MHD runs with ZEUS (From |5|) 

the SPH code to demonstrate that our results are truly independent of the details 
of the viscous dissipation, and so that our lack of an explicit viscosity does not 
affect our results. 

The kinetic energy decay curves for the four resolution studies are shown 
in Fig. ^ For each of our runs we performed a least-squares fit to the slope of 
the power-law portion of the kinetic energy decay curves. These results appear 
converged at the 5-10% level; it is very reassuring that the different numerical 
methods converge to the same result for the unmagnetized case. 

We find that highly compressible, isothermal turbulence decays close to lin- 
early in time, with rj — 0.98. Adding magnetic fields decreases the decay rate 
only slightly to -q ^^0. 85-0. 9, with very slight dependence on the field strength or 
adiabatic index. Similar results have been reported by Stone, Ostriker, & Gam- 
mie and Padoan & Nordlund |^ for compressible MHD turbulence, and by 
Biskamp & Miiller |^ for incompressible MHD turbulence. 

The clear astrophysical implication of these models is that even strong mag- 
netic fields, with the field in equipartition with the kinetic energy, cannot prevent 
the decay of turbulent motions on dynamical timescales. If molecular clouds live 
for longer than their dynamical time of roughly a megayear, as is generally be- 
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lieved even by those arguing for lifetimes under 10 Myr |^,^ rather than the 
more classical 30 Myr ||2^, then the significant kinetic energy observed in their 
gas must be supplied more or less continuously. 

3 Driven Turbulence 

To compute the energy dissipation from uniformly driven turbulence we initialize 
the turbulent flow with a narrow band of k values, using a top-hat function with 
roughly the same behavior as the steeply peaked curve used by Stone et al. Q in 
most of their models. We set the velocity field up as described above To drive 
the turbulence, we then normalize this fixed pattern to produce a set of pertur- 
bations 5i^{x, y, z), and at every time step add a velocity field 5v{x, y, z) = ASu 
to the velocity with the amplitude A now chosen to maintain constant kinetic 
energy input rate -Ejn = AE/ At. For compressible flow with a time-dependent 
density distribution, maintaining a constant energy input rate requires solving 
a quadratic equation in the amplitude A at each time step. For a grid with N 
zones on a side, each of volume ZiV^, the equation for A is 

1 ^ 

AE= -AV ^ pijkASvijk ■ (v^jk + ASu.jk)- (2) 

i,j,k—l 

We take the larger root of this equation to get the value of A. The resulting flow 
is shown in Fig. ^. 

We find that the best description of these compressible models comes by tak- 
ing a length scale C = Xd the driving wavelength, and a velocity scale V = Wrms 
the rms velocity, rather than any of the other options available. Happily, these 
are just the length and velocity scales that would be expected from the the- 
ory of incompressible turbulence. Figure ^ shows equilibrium energy dissipation 
rates for all the models described by Mac Low |5^, compared to the quantity 
^d^rms ^ '^rms/'^d- ^ to the hydrodyuamic models gives a relation with slope 
1.02. Let us define a dimensionalized wavenumber k = (27: / L)kd — 2'K/\d. A 
very good approximation is then the linear relation 

Ekin ^ -rjytnkv^^^, (3) 

with rjy = 0.21/7r, where the assumption is made that in equilibrium i?kin = £'in- 
The dependence on the mass of the cube m comes strictly from dimensional 
arguments, as all of the runs shown have the same mass m — pqL^ — 8. The 
referee of this review points out that rjykd — 0.42A^'^, remarkably close to the 
value of l/{2Xd) that would be expected by direct appHcation of the theory of 
incompressible turbulence. This implies that the density variations average out 
almost independently of velocity variations. 

The MHD models that fit equations most closely are the strong field 
cases, with vaJcs = 10. The weak field cases appear to follow a relation similar 
to equation (JSj), but with values of ?]„ up to a factor of two higher, as shown 
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Fig. 3. Models showing the appearance of turbulence with different characteristic size 
scales and driving energy input. Log of density on slices through hydrodynamic models 
from at a resolution of 128"^ grid points. The value of "drive" given in the figure is 
the driving energy input rate Ein for that model, while k = L/Xd is the wavenumber 
of the driving pattern, and the size of the cube L — 2 for all runs 



in Figure |^. Without further computation, it remains unclear how much of the 
variation seen among the models with the same /? is due to random fluctuations 
or the remaining lack of numerical convergence, and how much is real. The 
higher dissipation seen in the high-/?, weak-field cases could be explained by 
noting that weak fields will be more strongly infiucnced by the fiow, generating 
more dissipative MHD waves. 
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Fig. 4. Energy dissipation rate for hydro and MHD models compared to fcdWrms where 
kd is the dimensionless driving wavenumber, and the size of the cube L = 2 for all 
runs. The lines have slope of unity, and are ht to the hydro models, yielding the value 
for the dissipation coefficient rf^ = 0.21/7r (see equation Q VHydrodynamical models 
are indicated by squares, MHD models by triangles (From |5J) 



We can use these results to discuss whether decaying turbulence can delay 
gravitational collapse. This can formally be examined by determining whether 
the ratio 

r = td/ts > 1, (4) 

where the turbulent decay time td = -EkinZ-E'kin, and the free-fall time ts for the 
gas is given by equation (|l|) . Because td depends not only on the strength of the 
turbulence, but also on the driving wavelength, the value of r also depends on 
the ratio 

K = Xd/Xj, (5) 

where the Jeans wavelength Xj = Cs^/ir/Gpo. Numerical models that I disucss 
in the next section show that turbulence cannot support the gas against collapse 
at wavelengths longer than the driving wavelength |^^^, so that k < 1 is 
necessary for the turbulence to fully support against collapse. 
Substituting for the values in equation (^), we can write 

-Ekin Cs /32 

We can now use equation (^) for ii^kin, and, somewhat less accurately, take Skin ~ 
noting that this introduces no more than a 20-30% error. Substituting 
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Fig. 5. Dependence of the dissipation coefficient rjv — ii'kin/(fc'Ti'U?ms) on the plasma /3, 
the ratio of thermal to magnetic pressure. The dashed line shows the value derived from 
unmagnetized models. Models with weaker fields appear to have as much as a factor 
two higher dissipation rate, while strong field models approach the unmagnetized rate 
(From H) 

and using the definition of k given in equation (j^), I find that the dissipation 
time scaled in units of the free fall time is 

t(k) = — ~ 3.9 , (7) 

where M„„s = Uims/cs is the rms Mach number of the turbulence. In molecular 
clouds, Mrms IS typically observed to be of order 10 or higher, while k < 1 is 
required for support, so turbulence will decay long before the cloud collapses 
and not markedly influence its collapse. 



4 Self-Gravitating Turbulence 

Now that it is clear that driven turbulence must be present in molecular clouds, 
we can ask whether it alone is sufficient to support clouds against gravitational 
collapse. 

The virial theorem provides a first way of examining this question. In equi- 
librium the total kinetic energy in the system adds up to half its potential en- 
ergy, E'kin + l/2£'pot = 0. If i?kin + l/2i?pot < the system collapses, while 
£'kin + l/2£'pot > implies expansion. In turbulent clouds, the total kinetic 
energy includes not only the internal energy but also the contributions from 
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turbulent gas motions. If this is taken into account, simple energy considera- 
tions can already provide a qualitative description of the collapse behavior of 
turbulent self-gravitating media Q . 

A more thorough investigation, however, requires a linear stability analysis. 
For the case of an isothermal, infinite, homogeneous, self-gravitating medium at 
rest (i.e. without turbulent motions) Jeans |Q derived a relation between the 
oscillation frequency w and the wave number k of small perturbations, 

-c,^fc2 +47rGpo = , (8) 

where Cg is the isothermal sound speed, G the gravitational constant, and po the 
initial mass density. Note that the derivation includes the ad hoc assumption that 
the linearized version of the Poisson equation describes only the relation between 
the perturbed potential and the perturbed density, neglecting the potential of 
the homogeneous solution. This is the so-called 'Jeans swindle'. Perturbations 
are unstable against gravitational contraction if their wave number is below a 
critical value, the Jeans wave number fcj, i.e. if 

k^<k]^'-^, (9) 

or equivalently if the wave length of the perturbation exceeds a critical size given 
by A J = 2'KkJ^. This directly translates into a mass limit. All perturbations with 
masses exceeding the Jeans mass. 



PoX' = {^y^%o'^'cl (10) 



will collapse under their own weight. 

These and subsequent analytical approaches (reviewed in |^) make a strong 
assumption that substantially limits their reliability, namely that the equilib- 
rium state is homogeneous, with constant density po- However, observations 
clearly show that molecular clouds are extremely non-uniform. One way to 
achieve progress and circumvent the restrictions of a purely analytical approach 
is to perform numerical simulations. Bonazzola et al. |24| , for example, used 
low resolution (32 x 32 collocation points) calculations with a 2-dimensional 
spectral code to support their analytical results. Also restricted to two dimen- 
sions were the hydrodynamical studies by Passot et al. Leorat, Passot & 
Pouquet [Q, Vazquez-Semadeni et al. and Ballesteros-Paredes, Vazquez- 
Semadeni & Scalo [||, although performed with far higher resolution. Magnetic 
fields were introduced in two dimensions by Passot, Vazquez-Semadeni, & Pou- 
quet 1?^ , and extended to three dimensions with self-gravity (though at only 64'^ 
resolution) by Vazquez-Semadeni, Passot, & Pouquet |9^. A careful analysis of 
1-dimensional computations including both MHD and self-gravity was presented 
by Gammie & Ostriker |Q, who extended their work to 2.5 dimensions more 
recently | |6^ . 

We use ZEUS-3D and SPH to examine the gravitational stability of three- 
dimensional hydrodynamical turbulence at higher resolution than before, and 
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include magnetic fields using ZEUS-3D. The use of both Lagrangian and Eule- 
rian numerical methods to solve the equations of self-gravitating hydrodynamics 
in three dimensions (3D) allows us to attempt to bracket reality by taking ad- 
vantage of the strengths of each approach. This also gives us some protection 
against interpreting numerical artifacts as physical effects. 

SPH can resolve very high density contrasts because it increases the particle 
concentration, and thus the effective spatial resolution, in regions of high density, 
making it well suited for computing collapse problems as shown in Fig. ^. By the 




Fig. 6. SPH density cubes for model driven in the wavenumber interval 3 < fed < 4, 
shown (a) at the time when gravity is turned on, (b) when the first dense cores are 
formed and have accreted M« = 5% of the mass, (c) when the mass in dense cores is 
M, = 25%, and (d) when M« = 50%. Time is measured in units of the global system 
free-fall time scale Tff (From 
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same token, though, it resolves low-density regions poorly. Shock structures tend 
to be broadened by the averaging kernel in the absence of adaptive techniques. 
The correct numerical treatment of gravitational collapse requires the resolution 
of the local Jeans mass at every stage of the collapse In the current code, 
once an object with density beyond the resolution limit of the code has formed 
in the center of a collapsing gas clump it is replaced by a 'sink' particle p!7| . 
Adequately replacing high-density cores and keeping track of their further evo- 
lution in a consistent way prevents the time step from becoming prohibitively 
small. We are thus able to follow the collapse of a large number of cores until 
the overall gas reservoir becomes exhausted. 

ZEUS-3D, conversely, gives us equal resolution in all regions, and allows us 
to resolve shocks well everywhere, as well as allowing the inclusion of magnetic 
fields, as shown in Fig. 0. On the other hand, collapsing regions cannot be fol- 
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Fig. 7. Two-dimensional slices of 256'' models from driven with wavenumber k = 2 
hard enough that the mass in the box represents only 1/15 (Afj)turb, and with magnetic 
fields strong enough to give critical mass fractions as shown, corresponding to /9 = oo, 
0.9, and 0.04. Slices are taken at the location of the zone with the highest density at 
the time when 10% of the total mass has been accreted onto cores. The plot is centered 
on this zone. Arrows denote velocities in the plane. The length of the largest arrows 
corresponds to a velocity of u ~ 20cs. The density greyscale is given in the colorbar. 
(From S) 



lowed to scales less than one or two cells. We must again consider the resolution 
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required for gravitational collapse. For a grid-based simulation, the criterion 
given by Truelove et al. ||8^ holds. Equivalent to the SPH resolution criterion, 
the mass contained in one grid zone has to be rather smaller than the local Jeans 
mass throughout the computation. 

Applying this criterion strictly would limit our simulations to the very first 
stages of collapse, as we have not implemented anything like sink particles in 
ZEUS. We have therefore extended our models beyond the point of full resolution 
of the collapse, as we are primarily interested in the formation of collapsed 
regions, but not their subsequent evolution. Thus, in the ZEUS models, the fixed 
spatial resolution of the grid implies that strongly collapsed cores have a larger 
cross-section than appropriate for their mass. In encounters with shock fronts 
the probability for these cores to get destroyed or lose material is overestimated. 
Cores simulated with ZEUS are therefore more easily disrupted than they would 
be physically. SPH, on the other hand, underestimates the disruption probability, 
because sink particles cannot lose mass or dissolve again once they have formed. 
The physical result is thus bracketed by these two numerical methods. 

We use the same driving method as described above for both the SPH and 
ZEUS models. We define an effective turbulent Jeans mass (A/j)tuib by sub- 
stituting Cg — > Cg + 1/3 (u^) for the thermal sound speed Cg in equation ( p^ 
where we approximate the rms velocity of the flow (u^) by 2E]^inlM. The tur- 
bulent Jeans mass (Mj)turb must be compared to the total system mass M = 1 
in order to determine whether global stability is reached. 

We find that local collapse occurs even when the turbulent velocity field car- 
ries enough energy to counterbalance gravitational contraction on global scales, 
as shown in Fig. ^. This confirms the results of two-dimensional (2D) and low- 
resolution (64"^) 3D computations with and without magnetic fields by Vazquez- 
Semadeni et al. [ pO|. An example of local collapse in a globally supported cloud 
is given in Figure |6|. The presence of shocks in supersonic turbulence drasti- 
cally alters the result from analytic models of incompressible turbulence, as was 
first noted by Elmegreen and studied numerically by Vazquez- Semadeni et 



al. |90 . The density contrast in isothermal shocks scales quadratically with the 
Mach number, so the shocks driven by supersonic turbulence create density en- 
hancements with Sp oc Ai^, where M is the rms Mach number of the flow. In 
such fluctuations the local Jeans mass is decreased by a factor of M and therefore 
the likelihood for gravitational collapse increased. 

To test this explanation numerically, we designed a test case driven at short 
enough wave length and high enough power to support even fluctuations with 
dp cx A4^, and ran it with both codes, driven with a wave number k — 39 — 
40. Within 20 Tff this model shows no signs of collapse. All the other globally 
supported models with less extreme parameters that we computed did form 
dense cores during the course of their evolution, supporting our hypothesis that 
local collapse is caused by the density fluctuations resulting from supersonic 
turbulence. 

Magnetic fields might alter the dynamical state of a molecular cloud suf- 
ficiently to prevent gravitationally unstable regions from collapsing M. They 
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Fig. 8. Fraction of mass Af, in dense cores as function of time. All models are computed 
using SPH with sink particles replacing dense, collapsed cores. The different model 
parameters from j47j are indicated in the figure. The mass in the box is initially unity, 
so the solid curves are unsupported formally, while the others are supported. The figure 
shows how the efficiency of local collapse depends on the scale and strength of turbulent 
driving. Time is measured in units of the global system free-fall time scale rg, while 
the driving wavelength Ad = 27r/fc (From |M) 



have been hypothesized to support molecular clouds either magnetostatically or 
dynamically through MHD waves. Mouschovias & Spitzer derived an expres- 
sion for the critical mass-to-flux ratio in the center of a cloud for magnetostatic 
support. Assuming ideal MHD, a self-gravitating cloud of mass M permeated 
by a uniform flux (p is stable if the mass-to-flux ratio 



with Ccj depending on the geometry and the field and density distribution of 
the cloud. A cloud is termed subcritical if it is magnetostatically stable and 
supercritical if it is not. Mouschovias & Spitzer |]63| determined that c<j = 0.13 
for a spherical cloud. 

We include magnetic fields in our models of driven, self-gravitating turbu- 
lence to test their effectiveness in supporting against self-gravity. The MHD 
simulations start with a uniform magnetic fleld in the z-direction. We must con- 
sider the resolution required to accurately follow magnetized collapse. Numerical 
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diffusion can reduce the support provided by a static or dynamic magnetic field 
against gravitational collapse. Increasing the numerical resolution decreases the 
scale at which numerical diffusion acts. For strong, subcritical fields, the reso- 
lution should ensure that numerical diffusion remains unimportant even for the 
dense, shocked regions. We did a suite of models varying the sound speed and 
the mass in the cube while holding the magnetic field strength constant, thus 
varying M/Mcr and the number of zones in a Jeans length Aj. From these mod- 
els we conclude that for a self-gravitating magnetostatic sheet to be well resolved 
with the algorithm under study, its Jeans length must exceed four zones |^ . 

Collapse occurs in both unmagnetized and magnetized cases at all resolutions, 
as shown in Figure ^. However, increasing the resolution makes itself felt in 
different ways in hydrodynamical and MHD models. In the hydrodynamical case, 
higher resolution results in thinner shocks and thus higher peak densities. These 
higher density peaks form cores with deeper potential wells that accrete more 
mass and are more stable against disruption. If we increase the resolution in the 
MHD models, on the other hand, we can better follow short wavelength MHD 
waves, which appear to be able to delay collapse, although not to prevent it. 
This result extends to models with 512'^ zones [Q. 

5 Characterization of Turbulence 
5.1 Wavelet Transforms 

To obtain clues to the true physical nature of interstellar turbulence, character- 
istic scales and any inherent scaling laws have to be measured and modelled. A 
major problem with characterizing both the observations and the models is to 
determine what scaling behaviour, if any, is present in complex turbulent struc- 
tures. Both the velocity and density fields need to be considered, but only the 
radial velocity and column densities can be observed. 

One measure useful for characterizing structure and scaling in observed maps 
of molecular clouds is the zi-variance, a\ , an averaged wavelet transform method 
for measuring the amount of structure at different scales introduced by Stutzki et 
al. The Z\- variance spectrum clearly shows characteristic scales and scaling 
relations, and its logarithmic slope can be analytically related to the spectral 
index of the corresponding Fourier power spectrum (see Fig. [l^). In order to 
understand the physical significance of the characterization of the observational 
maps by Z\-variance spectra, we apply the same analysis to observations of the 
Polaris Flare at three different scales |]30| , |l9p^ ] , as analyzed by and to sim- 

ulated observations of gas distributions resulting from our MHD models (5^,^ , 
as well as to the actual density distributions from those models. 

In Fig. we compare the power spectrum and zi-variance for the density 
distibution from three simulations with different driving scales. Each shows a 
different characteristic scale formed by the driving, visible as a turn-over at 
large lags in the Z\-variances and at small wavenumbers in the power spectrum, 
respectively. At smaller lags and higher wavenumbers power laws can be seen 
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Fig. 9. Comparison of the mass accretion behaviour for runs driven at wavenumber 
fc = 1 — 2 with varying resolution. Pure hydro runs are shown in the upper panel, and 
MHD runs in the lower panel, with resolutions 64^ (dotted), 128^ (dashed), and 256^ 
(solid). M* denotes the sum of masses found in all cores determined by the modified 
clumpfind algorithm Note that these cores are subject to destruction as they 

can only collapse to the grid scale. Times are given in units of free-fall time. Although 
the collapse rate varies, we get collapse in all cases (From [53) 



in both cases for the models with longer wavelength driving, with their slopes 
related by the analytic relation mentioned above. A steep drop-off follows at the 
smallest scales indicating the resolution limit of the simulation. The model with 
the shortest wavelength (k^ — 8) driving shows no power law because the driving 
scale of 16 zones is only marginally larger than the dissipation scale of about 10 
zones. 

The Z\-variance analysis of astronomical maps was extensively discussed by 
Bensch et al. pO]. In all the observations they analyzed, the total cloud size was 
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Fig. 10. Comparison of 3D Zi-variance spectra (above) to Fourier power spectra (below) 
for 128"^ models of turbulence driven at wavenumbers of fed = 2 (solid), fed = 4 (dotted), 
and kd = 8 (dashed), demonstrating that the analytic relation between power spectrum 
slope and Zi-variance spectrum slope also holds for the local behavior of structures not 
showing a straight power-law (From j56|) 

the only characteristic scale detected by means of the Z\-variance. Below that 
size they found a self-similar scaling behaviour reflected by a power law with 
index a =0.5-1.3 corresponding to a Fourier power spectral index C =2.5-3.3. 

In Fig. in we show how numerical resolution, or equivalently the dissipation 
scale, influences the Z\-variance spectrum that we find from our simulations. We 
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test the influence of the numerical resolution on the structure by comparing a 
simple hydrodynamic problem of decaying turbulence computed at resolutions 
from 64'^ to 256'^, with an initial rms Mach number M = 5. In contrast to the 



0.1000 



0.0100 

1 

b 

0.0010 



1 


1 ' ' ' ' 












128^ 




1 


256^ 

1 . . . . 





0.01 0.10 1.00 

As 



Fig. 11. z4- variance of the density cubes (projected to two-dimensions) for a model of 
decaying hydrodynamic turbulence at a time t/ts = 0.1, where ts is the sound-crossing 
time, computed using numerical resolutions of 64'', 128'', and 256^ zones (From ]56|) 



results from which showed little dependence of the energy dissipation rate on 
the numerical resolution, we find here clear differences in the scaling behaviour 
of the turbulent structures. At small scales we find a very similar decay in the 
relative structure variations up to scales of about 10 times the pixel size (0.03, 
0.06, and 0.1 for the resolutions 256'^, 128^, and 64^, respectively) in all three 
models. This constant length range starting from the pixel scale clearly identifies 
this decay as being due to the numerical viscosity acting at the smallest available 
size scales. 

Similar behaviour can be observed at the largest lags where structure vari- 
ations decay for all three simulations at a length scale of roughly a quarter 
the cube size. This structure reflects the original driving of the turbulence with 
wavenumbers 1 < kd < 8, which produces structure on the corresponding length 
scales. Structures larger than about half the cube size are suppressed by the 
periodicity of the simulation cubes. 

We can now attempt to interpret the observations using the behavior seen in 
the model Z\-variance spectra, and equivalent behavior seen in velocity centroid 
Z\- variance spectra of the models, as described by Ossenkopf & Mac Low 
Fig. |l2| shows the square root of the Z\-variance for the Polaris Flare velocity 
centroid maps. An upturn in power at the smallest lags in each map produced by 
the observational noise adding power at small scales was subtracted to prepare 
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centroid A-variance 




Fig. 12. The square root of the Z\-variance of the smoothed Polaris Flare velocity 
centroid maps. The size of the smoothing filter was adapted to remove only the average 
noise level. Spectra are shown from three nested maps obtained with the 1.2 m CfA 
telescope with half power beam width (HPBW) 8.7', the 3 m KOSMA with HPBW 
2.2', and the 30 m IRAM with HPBW 0.35'. At a distance of 150 pc, 1000" = 0.73 pc. 
(From H) 



these spectra [20 . The velocity structure in the Polaris Flare maps shows power- 
law behaviour from scales of about 0.05 to several parsecs. Below that range 
the slope of the Zi- variance spectrum of the velocity centroids steepens, perhaps 
because of physical dissipation through ambipolar diffusion or other mechanisms, 
while above that range it flattens because of the lack of emission on larger scales. 



Bensch et al. [20 demonstrated similar scaling behaviour in the intensity maps, 
reflecting the density structure. 

Comparison between the observations of the Polaris Flare and the turbulence 
simulations constrains the mechanisms driving the turbulence in this cloud. Any 
mechanism that drives at an intermediate length scale, such as jets from embed- 
ded protostars, should produce a peak in the Z\-variance at that scale, which is 
not observed. The approximately self-similar, power-law behaviour seen in the 
observations is best reproduced by models where the energy is injected at large 
scales and dissipated at small scales. Driving by interactions with superbubbles 
and field supernova remnants p^ , p4[ would provide such a driving mechanism. 
The Polaris Flare molecular cloud lies in the wall of a large cylindrical struc- 
ture representing one of the nearest H I supershells, the North Celestial Polar 
Loop @, adding additional support to this proposal. 

The dominant physical mechanism for dissipation in molecular clouds was 
first shown by Zweibel & Josafatsson lol] to be ambipolar diffusion. Klessen et 
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al. showed that the length scale on which ambipolar diffusion will become 
important can be found by examining the ambipolar diffusion Reynolds number 

Ra = MaLv^/va (12) 

defined by Balsara ||l^ and Zweibel & Brandenburg where L and M.a are 
the characteristic length and Alfven Mach number, Vni = 7Pi is the rate at 
which each neutral is hit by ions, and v\ = B^/Anp approximates the effective 
Alfven speed in a mostly neutral region with total mass density p = Pi -\- Pn and 
magnetic field strength B. The coupling constant depends on the cross-section 
for ion-neutral interaction, and for typical molecular cloud conditions has a value 
of 7 « 9.2 X 1013 cm3 s'^ g-^ (e.g. Q). 

Setting the ambipolar diffusion Reynolds number Ra — 1 yields a diffusion 
length scale of 

Ld = VA/MAl^ni (13) 

/ B \ / 10-6 \ / 103 cm-3\'/' 

with the ionization fraction x — pi/pn and the neutral number density n„ = 
pnl with p = 2.36m//. If the ionization level in the Polaris Flare is low enough 
and the field is high enough, this length scale of order 0.05 pc would be directly 
resolved in the IRAM observations. We cannot yet unambiguously say whether 
the steepening observed at the smallest scales in the velocity spectra is due to 
beam smearing at the observational resolution limit or to a detection of the 
actual dissipation scale, similar to the downturn at the dissipation scale in the 
numerical models. If better observations do continue to show such a downturn 
in the future, that will be an indication of the dissipation scale. 



5.2 Clump Characterization 

Using 2D numerical simulations, Ballcsteros-Paredes et al. ||] showed that ob- 
served clumps frequently come from the superposition of several physically dis- 
connected regions in the line of sight at the same radial velocity (see also p7|) 
but not necessarily at the same position or three-dimensional velocity. This had 
earlier been suggested by workers including Burton ||2^, Issa et al. Falgar- 
one, Puget, & Perault and Adler & Roberts ||l|. 

We used the models described above to examine this effect in 3D using sim- 
ulated observations. We calculate the radiative transfer through the density and 
velocity fields given by the computations, assuming local thermodynamic equi- 
librium (LTE) for the population of the molecular energy levels, which is a 
sufficiently good approximation to study qualitatively the effects of projection 
and superposition (velocity crowding) of structure in the cloud. In order to de- 
termine differences when observing the same region with different tracers, we 
also set minimum density thresholds below which the molecules might be under- 
excited. 
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Velocity crowding contributes to the generation of clumps in the observa- 
tional space, so observed clumps frequently contain emission from physically 
separated regions (8|j7|,|^,|9|. We demonstrate this effect using a typical MHD 
simulation with intermediate values of driving strength (L = 1) and wavenumber 
(k = 4). 

Figure |l^ shows physical and observed maps for ^■^CO(l-O) and (2-1). We 
see that clumps in real space (letters A, B, and C in panel y — x{l — 0)) do not 
necessarily have a counterpart in observational space. Clumps in observational 
space, on the other hand, (letter D in panel Vy—x{l — 0)) do not necessarily come 
from isolated regions in real space, but have contributions from many different 
regions along the same line of sight. In Figure |l^ we plot dotted lines that show 
the places where the emission of a physical clump lies in observational space, and 
where the emission from an observational clump is generated in physical space. 
For reference, we use the same lines in the ^'^CO(l-O) map as in the ^'^C0(2-l) 
map. 

If clumps in the observational space are the result of the contribution of 
multiple regions in the physical space it is of primary importance to understand 
whether relationships reported for observed clumps in molecular clouds are also 
valid for the actual physical clumps. 

Larson |Q studied the dependence with size of the mean density, velocity 
dispersion, and mass spectrum of the clouds in a sample of observational data 
taken from the literature. He found 

p cx i?", (15) 

Sv (X i?", (16) 
dN 



dlogM 



oc M'^. (17) 



The most commonly quoted values in the literature are a~ —1, m ^ 1/2, and 
7 ~ —0.5. However, there is some discrepancy in the values reported p^ , |5T[ . 

Larson himself mentioned that the relationships he found might be due 
to the limitations of the observations. Kegel first demonstrated that the ob- 
served mean density-size relationship could be due to observational effects, and 
that the observed and physical values of properties such as radius or volumetric 
density might be quite different. Scalo |Q showed that CO and extinction satu- 
rate at roughly the same column density, which forms the upper envelope of the 
mean density-size relationship. Vazquez-Semadeni et al. |Q reported the lack of 
a mean density-size relationship, confirming numerically the analysis by Kegel 
(01, see also the discussions in |Q and |Q) in the sense that there are clouds 
with small sizes and low column density that will be undetected in observational 
surveys. 

In Figure |lj we use the computational model described above to compare the 
density-size relationships in physical and observational spaces. We define clumps 
as a connected set of points below a local maximum following the intensity only 
downwards until the threshold is reached, a scheme implemented in the code 
called clumpf ind |92 . We note three points about Figure First, there is no 
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Fig. 13. Physical and velocity space maps for ^^CO(l-O) and (2-1) in a run with E = 1, 
fc = 4 and B — 0.1. Clumps in physical space (A, B, C) do not necessarily correspond 
to clumps in the observed velocity space. Observed clumps in velocity space (D) are 
not necessarily formed by emission from a single region in physical space. From M 

relation between mean density and size (eq. p^), confirming the 2D results of 
Vazquez- Semadeni et al. [^^. Second, there is a minimum density below which 
there are no clumps identified. This minimum is just given by the density thresh- 
old we used in clumpf ind. Third, even though the simulations exhibit a large 
dynamical range in density (/Omax/pmin ~ 3.5 x 10"*), the dynamical range in the 
mean density-size relationship is small, because in constructing such a plot we 
choose the clumps around the local density maxima. 

In the middle and lower panels of Fig. ^ we show the mean density-size 
relationship for clumps in simulated observational maps of the same run, inte- 
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Fig. 14. Mean density-size relationship for physical clumps in physical (PPP, upper 
panel); and simulated observational clumps in two different density tracers in obser- 
vational coordinates (PPV, middle and lower panels). The dotted line has a slope of 
a = —1. In physical space we find no correlation, verifying the results of Vazquez- 
Semadeni et al. [Q, but nevertheless the simulated observations show such a correla- 
tion, as found by Larson and many others. The selection of two different density 
tracers was chosen to show that the apparent correlation does not depend on the se- 
lection of the density threshold. (From |7|) 



grated along the x axis. The middle panel is the one obtained by using CS(l-O), 
and the lower panel is the relationship obtained by using ^'^CO(l-O). The ob- 
served clumps do exhibit approximately the relationship given by Larson |Q, 
despite the lack of correlation exhibited by the physical clumps in this model. 
We conclude from this demonstration that the observed density-size relationship 
(eq. H^) is an observational artifact. 
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Kegel suggested one mechanism that appears able to explain our results: 
only clouds with emission exceeding some intensity threshold will be detected, 
effectively setting a column density cutoff, rather than the physical density cutoff 
imposed in the physical density-size relationship. (This threshold is determined 
by the noise level in observations, or by some average intensity in the simulated 
observations). A constant column density cutoff produces a cutoff with slope —1 
in the mean density-size plane (middle and lower panels), just as the constant 
physical density cutoff in physical space produces a flat cutoff. 

6 Supernova-Driven Turbulence 

Both support against gravity and maintenance of observed motions appear to 
depend on continued driving of the turbulence, as I have described. What then 
is the energy source for this driving? 

Motions coming from gravitational collapse have often been suggested, but 
fail due to the quick decay of the turbulence as described above. If the turbulence 
decays in less than a free-fall time, then it cannot delay collapse for substantially 
longer than a free-fall time p6| . 

Protostellar jets and outflows are another popular suspect for the energy 
source of the observed turbulence. They are indeed quite energetic, but they 
deposit most of their energy into low density gas, as is shown by the observa- 
tion of multi-parsec long jets extending completely out of molecular clouds 
Furthermore, as described in the previous section, the observed motions show 
increasing power on scales all the way up to and perhaps beyond the largest 
scale of molecular cloud complexes It is hard to see how such large scales 
could be driven by protostars embedded in the clouds. 

Another energy source that has long been considered is shear from galactic 
rotation. Work by Sellwood & Balbus has shown that magnetorotational 
instabilities |^,^ could couple the large-scale motions to small scales efficiently. 
For parameters appropriate to the far outer H i disk of the Milky Way, they 
derive a resulting velocity dispersion of 6 km s~^, close to that observed. This 
instability may provide a base value for the velocity dispersion below which no 
galaxy will fall. If that is sufficient to prevent collapse, little or no star formation 
will occur, producing something like a low surface brightness galaxy with large 
amounts of H i and few stars. 

In active star-forming galaxies, however, clustered and field supernova ex- 
plosions appear likely to dominate the driving, raising the velocity dispersion to 
the 10-15 km s^^ observed in star-forming portions of galaxies (see work cited 
in ||5^ for example). These explosions will be predominantly from B stars no 
longer associated with their parent gas, as they are far more numerous than 
more massive O stars, have explosions that are just as powerful, and live long 
enough (up to 50 Myr) to either drift away from or ionize their parent clouds. 
This provides a large-scale self-regulation mechanism for star formation in disks 
with sufficient gas density to collapse despite the velocity dispersion produced by 
the magnetorotational instability. As star formation increases in such galaxies. 
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the number of OB stars increases, ultimately increasing the supernova rate and 
thus the velocity dispersion, which will restrain further star formation. 

Supernova driving not only determines the velocity dispersion, but may actu- 
ally form molecular clouds by sweeping gas up in a turbulent flow. Clouds that 
are turbulently supported will experience inefficient, low-rate star formation, 
while clouds that are too massive to be supported will collapse undergoing 
efficient star formation to form OB associations or even starburst knots. 

We study supernova driving in a magnetized medium numerically, using the 
RIEMANN framework for computational astrophysics, which is based on higher- 
order Godunov schemes for MHD [ ^6l|Tll , [l^ , and incorporates schemes for pres- 
sure positivity and divergence- free magnetic fields (The framework 
also includes parallelized, divergence-conserving, MHD adaptive mesh refine- 
ment [p^ljl^ , though no results using that capability are shown here.) In the 
models presented here, we solve the ideal MHD equations including both radia- 
tive cooling and pervasive heating in a (200 pc)'^ periodic computational box, 
using a grid of 128^ cells. We start the simulations with a uniform density of 
2.3 X 10"^"* g cm~'^, threaded by a uniform magnetic field in the z-direction 
with strength 5.8 fiG, a factor of roughly two stronger than that observed in the 
Milky Way disk. This very strong field maximizes the effects of magnetization 
on the turbulence. 

For the cooling, we use a tabulated version of the radiative cooling curve 
shown in Figure 1 of MacDonald and Bailey |5^, which is based on the work 
of Raymond, Cox & Smith |^ and Shapiro and Moore [Q. (It falls smoothly 
from temperatures of order lO"" K to 10^ K, not incorporating a sharp cutoff at 
10* K due to the turnoff of Lya cooling.) In order to prevent the gas from cooling 
below zero, we set the lower temperature cutoff for the cooling at 100 K. We 
also include a diffuse heating term to represent processes such as photoelectric 
heating by starlight, which we set constant in both space and time. We set the 
heating level such that the initial equilibrium temperature determined by heating 
and cooling balance is 3000 K. 

We explode SNe at a rate of one every 0.1 Myr in our box, twelve times 
higher than our present Galactic rate, corresponding to a mild starburst like 
M82. The SNe are permitted to explode at random positions. Each SN explosion 
dumps 10^^ erg thermal energy into a sphere with radius 5 pc. The evolution 
of the system is determined by the energy input from SN explosions and diffuse 
heating and the energy lost by radiative cooling. An example of the resulting 
turbulence is shown in Figure 

The first theories of the multi-phase ISM, such as Field, Goldsmith, & Habing j 
postulated an isobaric medium. Since then, multi-phase models have commonly 
been interpreted as being isobaric, although McKee & Ostriker and Wolfire 
et al. ||9^ actually assume only local pressure equilibrium, not global, and McKee 
& Ostriker [s^ considered the distribution of pressures. In typical multi-phase 
models, the heating and cooling rates of the gas have different dependences on 
the temperature and density, so that the balance between heating and cooling 
determines allowed temperatures and densities for any particular pressure. This 
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Fig. 15. Two-dimensional slices through the three-dimensional MHD model M2, par- 
allel to the magnetic field at a time of 6.6 Myr, showing density (upper left), thermal 
pressure (upper right), temperature (lower left), and magnetic pressure (lower right). 
Color bars indicate the scale of each quantity. (From |b4]) 



balance can be shown graphically in a phase diagram, showing, for example, the 
allowed densities for any pressure (|Q: for a modern example, see Fig. 3(a) of 

In Figure |16|, the thermal-equilibrium curve for the heating and cooling mech- 
anisms included is shown as a black line. Only a single phase is predicted at high 
densities as our cooling curve did not include the physically-expected unstable 
region at temperatures of order 10^ K [0. Thus, if our model produced an 
isobaric medium, it would be expected to have a single low-temperature phase 
in uniform density given by the point at which the thermal-equilibrium curve 
crosses that pressure level. (Effectively, we would have the hotter two of the 
three phases proposed by McKee & Ostriker ||5^.) 




Fig. 16. Scatter plot of pressure vs. density at t = 6.6 Myrs in the MHD simulation 
M2, showing 32"^ points sampled at intervals of four points in each direction. Note that 
for each density a wide variation in pressure is seen. Cool gas with logT < 3.9 is shown 
in blue, warm gas with 3.9 < logT < 5.5 in green, and hot gas with logT > 5.5 in red. 
The thermal equilibrium curve for the cooling and heating functions in this simulation 
is overlaid as a black line. (The line of points at the very bottom right corresponds to 
an absolute cutoff in the cooling at 100 K that was enforced on the temperature in this 
model.) Note that our cooling curve may artificially prevent much low-temperature gas 
from forming in this model. From H] 



The scattered points in Figure |16| show the actual density and pressure of 
individual zones in the model. Many zones at low temperature do lie on the 
thermal equilibrium curve, but scattered all up and down it at many different 
pressures and densities. Furthermore, a substantial fraction of the gas has not 
had time to reach thermal equilibrium at all after dynamical compression. It 
appears that pressures are determined dynamically, and the gas then tries to 
adjust its density and temperature to reach thermal equilibrium at that pressure. 
Most gas will land on the thermal equilibrium curve when dynamical times are 
long compared to heating and cooling times. 

Even if we include the proper physics to allow multiple phases, the behav- 
ior observed in our model will lead to all points within the range of pressures 
available along the thermal equilibrium line being occupied, rather than the ap- 
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pearance of discrete phases (also see the chapter by Vazquez-Semadeni et al. in 
this volume). Unstable regions along the thermal equilibrium curve and off 
it will also be populated, as observed by Heiles |3^, but not as densely, as gas 
will indeed attempt to heat or cool to a stable thermal equilibrium at its current 
pressure. In particular, cold high pressure regions can be formed dynamically, 
without the influence of self-gravity, perhaps giving a method for forming molec- 
ular clouds with observed properties that are not in hydrostatic equilibrium. 

The relation between magnetic and thermal pressure is shown in Figure 
In this Figure, the relative strength of thermal and magnetic pressure is shown 
at one time for the magnetized simulation. The scattering of regions at very low 
thermal pressure all have substantial magnetic pressures, demonstrating that 
magnetically supported regions can occur. However, their relative importance 
is rather low, as shown by the small number of points in that regime. Hot gas 
can be seen, on the other hand, to be dominated by thermal pressure, with low 
magnetic pressures. 

As we have shown, the supernova-driven models have broad ranges of pres- 
sures. We can quantify this by examining the pressure probability density func- 
tion (PDF), as shown in Figure [T8|. In both cases, these show roughly log-normal 
pressure PDFs, very unlike the power-law distributions predicted by the analytic 
theory derived by McKee & Ostriker ||5^ . The observed distributions rather more 
resemble the log-normal density distributions suggested by Passot & Vazquez- 
Semadeni 1?^. Mac Low et al. show that the cool gas can actually be modeled 
quite successfully with this heuristic theory. 

7 Conclusions 

• Even relatively strong magnetic fields, with the field in equipartition with the 
kinetic energy, cannot prevent the decay of turbulent motions on dynamical 
timescales far shorter than the observed lifetimes of molecular clouds. The 
significant kinetic energy observed in molecular cloud gas must be supplied 
more or less continuously. 

• Supersonic turbulence strong enough to globally support a molecular cloud 
against collapse will usually cause local collapse. The turbulence establishes 
a complex network of interacting shocks. The local density enhancements in 
fluctuations created by converging shock flows can be large enough to become 
gravitationally unstable and collapse. The probability for this to happen, the 
efficiency of the process, and the rate of continuing accretion onto collapsed 
cores are strongly dependent on the driving wave length and on the rms 
velocity of the turbulent flow, and thus on the driving mechanism. 

• Interstellar clouds driven on large scales or without even global turbulent 
support very rapidly form stars in clusters. On the contrary, in gas that is 
supported by turbulence, local collapse occurs sporadically over a large time 
interval, forming isolated stars. The total star formation efficiency before the 
cloud dissolves due to stellar feedback or external shocks will probably be 
low. Thus, the strength and nature of the turbulence may be fully sufficient 
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Fig. 17. Scatter plot of magnetic vs. thermal pressure at i = 6.6 Myrs in the MHD 
simulation. We again plotted a subset of 32^ points sampled at intervals of four points 
in each direction. Note that regions of very low thermal pressure have substantial 
magnetic pressures. From M 



to explain the difference between the observed isolated and clustered modes 
of star formation. 

• Magnetorotational instabilities may provide a base value for the velocity 
dispersion below which no galaxy will fall. If that is sufficient to prevent 
collapse, little or no star formation will occur, producing something like a 
low surface brightness galaxy with large amounts of H i and few stars. In 
star-forming galaxies, however, clustered and field supernova explosions, pre- 
dominantly from B stars no longer associated with their parent gas, appear 
likely to dominate the driving, raising the velocity dispersion to some 10-15 
km s^^. 

• In a supernova-driven interstellar medium, we find a broad range of pres- 
sures with a log-normal distribution, and a substantial fraction of associated 
densities far from the thermal equilibrium values. This limits the predictive 
usefulness of phase diagrams based on thermal equilibrium, although thermal 
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Fig. 18. (6) Volume-weighted PDFs of pressure from the MHD models for cool gas 
with logT < 3.9 (blue), warm gas with 3.9 < logT < 5.5 (green), and hot gas with 
logT > 5.5 (red) at different times of 5.55 Myr (dashed), 6.06 Myr (dotted), and 
6.55 Myr (solid) in the 1.6 pc resolution model M2. This model has twelve times the 
galactic rate of SNe, which results in a very broad pressure distribution, (c) Mass 
distribution of pressure for the MHD model M2 for the full distribution (black), and 
for cool gas with logT < 3.9 (blue), warm gas with 3.9 < logT < 5.5 (green), and hot 
gas with logT > 5.5 (red) at a time of 6.55 Myrs. Most of the mass is found in cold 
gas, with a broad distribution around the peak pressure. From |m| 



equilibrium at the local pressure will still be the mildly favored state. Gas 
pressures appear to be determined dynamically, and each individual parcel of 
gas seeks local thermal equilibrium at the pressure imposed on it by the tur- 
bulent flow. Inferences that molecular clouds must be gravitationally bound 
because of their high observed confinement pressures are called into question 
by these results. Regions with densities approaching the overall densities of 
GMCs, and pressures an order of magnitude above the average interstellar 
pressure appear in our simulations even in the absence of self-gravity. 
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